#include <life/options.hpp>
#include <life/lifealg/backend.hpp>
#include <life/lifediscr/functionspace.hpp>
#include <life/lifediscr/region.hpp>
#include <life/lifepoly/im.hpp>
#include <life/lifefilters/gmsh.hpp>
#include <life/lifefilters/exporter.hpp>
#include <life/lifepoly/polynomialset.hpp>
#include <life/lifevf/vf.hpp>

#include <life/lifefilters/geotool.hpp>

/** use Life namespace */
using namespace Life;
using namespace Life::vf;


namespace MortarMethod
{

  inline
  po::options_description
  makeOptions()
  {
    po::options_description code1doptions("Code1D options");
    code1doptions.add_options()
      ("hsize", po::value<double>()->default_value( 0.1 ), "mesh size")
      ("hsize1", po::value<double>()->default_value( 0.1 ), "mesh size1")
      ("hsize2", po::value<double>()->default_value( 0.1 ), "mesh size2")
      ("nu", po::value<double>()->default_value( 0.5 ), "nu")
      ("alphaX", po::value<double>()->default_value( 0.01 ), "alphaX")
      ("alphaY", po::value<double>()->default_value( -2.0 ), "alphaY")
      ("bccoeff", Life::po::value<double>()->default_value( 100.0 ), "coeff for weak Dirichlet conditions")
      ("toler", Life::po::value<double>()->default_value( 1e-6 ), "error tol")
      ("xc", Life::po::value<double>()->default_value( 0.0 ), "param de recouvrement des maillages")
      ("penaldir", Life::po::value<double>()->default_value( 10 ),
       "penalisation parameter for the weak boundary Dirichlet formulation")
      ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh( 2.0 or 2.1)")
      ;
    return code1doptions.add( Life::life_options() );
  }

  inline
  AboutData
  makeAbout()
  {
    AboutData about( "MortarMethod2" ,
		     "MortarMethod2" ,
		     "0.2",
		     "Code MortarMethod",
		     Life::AboutData::License_GPL,
		     "Copyright (c) 2010 Universite Joseph Fourier");

    about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
    return about;

  }


  template<uint Dim,uint N,uint OrderGeo>
  class Localisation
    :
    public Application
  {
    typedef Application super;
  public:

    //! Polynomial order \f$P_2\f$
    //    static const uint16_type Order = 5;
    ///static const uint Order = N;

    typedef double value_type;

    //---------------------------------------------------------------------------------------//

    typedef Backend<value_type> backend_type;
    typedef boost::shared_ptr<backend_type> backend_ptrtype;

    typedef backend_type::sparse_matrix_type sparse_matrix_type;
    typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
    typedef backend_type::vector_type vector_type;
    typedef backend_type::vector_ptrtype vector_ptrtype;

    //---------------------------------------------------------------------------------------//

    typedef Simplex<Dim,OrderGeo,Dim> convex_type;
    typedef Mesh<convex_type> mesh_type;
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

    typedef Simplex<1,OrderGeo,2> convex_1d_type;
    typedef Mesh<convex_1d_type> mesh_1d_type;
    typedef boost::shared_ptr<mesh_1d_type> mesh_1d_ptrtype;

    //---------------------------------------------------------------------------------------//

    typedef bases<Lagrange<N,Scalar,PointSetFekete> > basis_u_type;
    typedef bases<Lagrange<0,Scalar> > basis_l_type;

    typedef FunctionSpace<mesh_type, basis_u_type> space_u_type;
    typedef boost::shared_ptr<space_u_type> space_u_ptrtype;

    typedef typename space_u_type::element_type element_u_type;

    typedef FunctionSpace<mesh_1d_type, basis_l_type> space_l_type;
    typedef boost::shared_ptr<space_l_type> space_l_ptrtype;

    typedef typename space_l_type::element_type element_l_type;

    //---------------------------------------------------------------------------------------//

    typedef Exporter<mesh_type> export_type;
    typedef boost::shared_ptr<export_type> export_ptrtype;

    /**
     * Constructor
     */
    Localisation( int argc, char** argv, AboutData const& ad,
		  po::options_description const& od )
      :
      super( argc, argv, ad, od ),
      //M_backend( backend_type::build( this->vm() ) ),
      //M_backend2( backend_type::build( this->vm() ) ),
      //M_backend( backend_type::build( BACKEND_TRILINOS ) ),
      meshSize( this->vm()["hsize"].template as<double>() ),
      meshSize1( this->vm()["hsize1"].template as<double>() ),
      meshSize2( this->vm()["hsize2"].template as<double>() ),
      exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) ),
      exporter1( Exporter<mesh_type>::New( this->vm(), "domain_1" ) ),
      exporter2( Exporter<mesh_type>::New( this->vm(), "domain_2" ) )
    {
    }

    mesh_ptrtype createMesh( std::string __name,std::string __str) ;
    mesh_1d_ptrtype createMesh1d( std::string __name,std::string __str) ;

    mesh_ptrtype createMesh( double meshSize );
    mesh_ptrtype createMesh1( double meshSize ,double x_couvr);
    mesh_ptrtype createMesh2( double meshSize, double x_couvr);

    mesh_1d_ptrtype createMeshInterface( double meshSize, double x_couvr);

    
    // run the application
    void runMultiDomain();

  private:

    // solve the system \f$D u = F\f$
    void solve( sparse_matrix_ptrtype const& D, element_u_type& U, vector_ptrtype const& F, bool is_sym );
    //void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
    void problem1(element_u_type& u , element_l_type& l  );
    void problem2(element_u_type& u , element_l_type& l  );

    // export results to ensight format (enabled by  --export cmd line options)
    void exportResults( element_u_type& U1,element_u_type& U2,double time );

    template<typename IteratorRange, typename QuadratureT, typename ExprT>
    void integrateLocal( IteratorRange const& itrange,
			 QuadratureT const& im,
			 ExprT const& expr,
			 vector_ptrtype& rhs );

  private:
    /*
    // linear algebra backend
    backend_ptrtype M_backend,M_backend2;
    sparse_matrix_ptrtype A1;
    vector_ptrtype B1;
    sparse_matrix_ptrtype A2;
    vector_ptrtype B2;
    */

    // mesh characteristic size
    double meshSize,meshSize1,meshSize2;

    // exporter factory
    export_ptrtype exporter;
    export_ptrtype exporter1;
    export_ptrtype exporter2;
  }; // Interpolate



  template<uint Dim,uint N,uint OrderGeo>
  typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
  Localisation<Dim,N,OrderGeo>::createMesh( std::string __name,std::string __str)
  {
    //typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
    //typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
    mesh_ptrtype mesh( new mesh_type );

    Gmsh gmsh;
    gmsh.setOrder(OrderGeo);

    std::string fname = gmsh.generate( __name, __str  );

    ImporterGmsh<mesh_type> import( fname );

    float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
    if (gmsh_version==2.0) import.setVersion( "2.0" );
    else if (gmsh_version==2.1) import.setVersion( "2.1" );

    mesh->accept( import );
    mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
    mesh->updateForUse();

    return mesh;

  }

  template<uint Dim,uint N,uint OrderGeo>
  typename Localisation<Dim,N,OrderGeo>::mesh_1d_ptrtype
  Localisation<Dim,N,OrderGeo>::createMesh1d( std::string __name,std::string __str)
  {
    //typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
    //typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
    mesh_1d_ptrtype mesh( new mesh_1d_type );

    Gmsh gmsh;
    gmsh.setOrder(OrderGeo);

    std::string fname = gmsh.generate( __name, __str  );

    ImporterGmsh<mesh_1d_type> import( fname );

    float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
    if (gmsh_version==2.0) import.setVersion( "2.0" );
    else if (gmsh_version==2.1) import.setVersion( "2.1" );

    mesh->accept( import );
    //mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
    mesh->updateForUse();

    return mesh;

  }


  template<uint Dim,uint Order,uint OrderGeo>
  typename Localisation<Dim,Order,OrderGeo>::mesh_ptrtype
  Localisation<Dim,Order,OrderGeo>::createMesh1( double meshSize,double x_couvr )
  {
    mesh_ptrtype mesh( new mesh_type );
    GeoTool::Node x1(-1,-1);
    GeoTool::Node x2(0,1);
    GeoTool::Rectangle R(meshSize,"Omega1",x1,x2);
    R.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
    R.setMarker(_type="line",_name="Interface",_marker2=true);
    R.setMarker(_type="surface",_name="Omega1",_marker1=true);
      
    return createMesh("Omega1",R.geoStr());

  }
  //----------------------------------------------------------------------------------------------------------//
  //----------------------------------------------------------------------------------------------------------//
  //----------------------------------------------------------------------------------------------------------//

  template<uint Dim,uint Order,uint OrderGeo>
  typename Localisation<Dim,Order,OrderGeo>::mesh_ptrtype
  Localisation<Dim,Order,OrderGeo>::createMesh2( double meshSize, double x_couvr )
  {
    mesh_ptrtype mesh( new mesh_type );
    GeoTool::Node x1(0,-1);
    GeoTool::Node x2(1,1);
    GeoTool::Rectangle R(meshSize,"Omega2",x1,x2);
    R.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker2=true);
    R.setMarker(_type="line",_name="Interface",_marker4=true);
    R.setMarker(_type="surface",_name="Omega2",_marker1=true);
      
    return createMesh("Omega2",R.geoStr());
  }

  template<uint Dim,uint Order,uint OrderGeo>
  typename Localisation<Dim,Order,OrderGeo>::mesh_1d_ptrtype
  Localisation<Dim,Order,OrderGeo>::createMeshInterface( double meshSize, double x_couvr )
  {

    mesh_1d_ptrtype mesh( new mesh_1d_type );

    std::ostringstream ostr;
    ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=" << meshSize << ";\n"
         << "Point(1) = { 0, -1, 0, h};\n"
         << "Point(2) = { 0, 1, 0, h};\n"
         << "Line(1) = {1, 2};\n"
         << "Physical Line(\"Interface\") = {1};\n";


      /*
	<< "Line(3) = {3, 4};\n"
	<< "Line(4) = {4, 1};\n"
         << "Line Loop(1) = {1, 2, 3 ,4 };\n"
         << "Plane Surface(2) = {1};\n"
         << "Physical Line(\"Force\") = {1};\n"
         << "Physical Line(\"Fixe\") = {2,4};\n"
         << "Physical Line(\"Libre\") = {3};\n"
         << "Physical Surface(\"Physique1\") = {2};\n";
      */
    return createMesh1d("OmegaInterface",ostr.str());
  }


  template<uint Dim,uint Order,uint OrderGeo>
  void
  Localisation<Dim,Order,OrderGeo>::exportResults( element_u_type& u1, element_u_type& u2, double time  )
  {
  
    exporter1->step( time )->setMesh( u1.functionSpace()->mesh() );
    exporter1->step( time )->add( "u1", u1 );
    //exporter1->step( time )->add( "lambda1", p1 );
    exporter1->save();

    exporter2->step( time )->setMesh( u2.functionSpace()->mesh() );
    exporter2->step( time )->add( "u2", u2 );
    //exporter2->step( time )->add( "lambda2", p2 );
    exporter2->save();

  } // Test_Schwarz::exportResults


  template<uint Dim,uint Order,uint OrderGeo>
  void
  Localisation<Dim,Order,OrderGeo>::runMultiDomain()
  {

    if ( this->vm().count( "help" ) )
      {
	std::cout << this->optionsDescription() << "\n";
	return;
      }

    using namespace Life::vf;  
 
    //-----------------------------------------------------------------------------------//

    //int MyOrder=N;

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			    % this->about().appName()
			    % convex_type::name()
			    % Order
			    % this->vm()["hsize"].template as<double>()
			    );

    //-----------------------------------------------------------------------------------//

    double xc = this->vm()["xc"].template as<double>();
    mesh_ptrtype mesh1 = createMesh1( meshSize1,xc );
    mesh_ptrtype mesh2 = createMesh2( meshSize2,xc );


    mesh_1d_ptrtype mesh1d = createMeshInterface(meshSize1,xc);

    space_u_ptrtype Xh1 = space_u_type::New( mesh1 );
    space_u_ptrtype Xh2 = space_u_type::New( mesh2 );

    space_l_ptrtype Lh = space_l_type::New( mesh1d );

    element_u_type u1( Xh1, "u1" );
    element_u_type u2( Xh2, "u2" );

    element_l_type l( Lh, "l" );

    //element_type U1old( Xh1, "U1old" );
    //element_type U2old( Xh2, "U2old" );


    //double __ppp =math::sqrt( integrate( elements(mesh1), _Q<2*Order>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
    //double __ppp2 =math::sqrt( integrate( elements(mesh2), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );

    //
  
    //-----------------------------------------------------------------------------------//

    /*A1=M_backend->newMatrix(Xh1,Xh1);
    B1=M_backend->newVector( Xh1 );

    A2= M_backend->newMatrix(Xh2,Xh2);
    B2= M_backend->newVector( Xh2 );
    */
    value_type nu=1;// this->vm()["nu"].template as<value_type>();
    AUTO( pi, M_PI );
    AUTO( Uex, cos(pi*Px()/2.)*cos(pi*Py()/2.) );
    AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
    AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );
    AUTO( f, -nu*lap_Uex);

    backend_ptrtype backend(backend_type::build( this->vm() ));

    sparse_matrix_ptrtype B1 = backend->newMatrix(Lh,Xh1);
    form2( Lh, Xh1, B1, _init=true ) = integrate( elements(mesh1d),idt(l)*id(u1) );

    sparse_matrix_ptrtype A1 = backend->newMatrix(Xh1,Xh1);
    vector_ptrtype F1 = backend->newVector(Xh1);

    form2( Xh1, Xh1, A1, _init=true ) = integrate( elements(mesh1),idt(u1)*id(u1) );
    form1( Xh1, F1, _init=true ) = integrate( elements(mesh1),f*id(u1) );

    form1( Xh1, F1) += integrate( markedfaces(mesh1, mesh1->markerName("Interface")),trans(grad_Uex)*N()*id(u1) );

    A1->close();
    F1->close();

    form2( Xh1, Xh1, A1 ) +=
      on( markedfaces(mesh1, mesh1->markerName("Bord")) , u1, F1, cst(0.) );

    //form2( Xh1, Xh1, A1 ) +=
    //  on( markedfaces(mesh1, mesh1->markerName("Interface")) , u1, F1, Uex );

    vector_ptrtype Am1F1= backend->newVector( Xh1 );
    backend->solve( A1, A1, Am1F1, F1, false );
    u1 = *Am1F1;

    exportResults( u1, u1, 0);

    //B2= M_backend->newVector( Xh2 );


    //-----------------------------------------------------------------------------------//




    //form2( Lh, Xh1, BB1, _init=true ) 
    //= integrate( markedfaces(mesh1,mesh1->markerName("Interface")),idt(l)*id(u1) );
    //form2( Lh, Lh, BBB1, _init=true ) = integrate( elements(mesh1d),-id(l)*idt(l)*id(u1) );


    //form1( Xh, B, _init=true );
      //integrate( elements(mesh), _Q<7>(),
      //f*id(v) );

    /*typedef typename _Q<3*Order>::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex>::type quad_type2;
    quad_type2 __quadra;
    this->integrateLocal( markedfaces( mesh1d, mesh1d->markerName( "Interface" )),
			  __quadra,//_Q<3*Order>(),
			  Px(),//grad_g*N() 
			  BB1 );
    */
    //-----------------------------------------------------------------------------------//

    AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );
  

#if 0
    double erreur=1;
    double err_tol=this->vm()["toler"].template as<double>();
    double err1,err2,norm1,norm2;
    int cpt=0;
    std::cout <<"\nola\n";    
    while (/*erreur>=err_tol&&*/ cpt <10) 
      {
   
	//U1old=U1;
	u1old=u1;
	p1old=p1;

	this->problem1(U1,U2);
	norm1 =math::sqrt( integrate( elements(mesh1), _Q<2*Order>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
	err1=math::sqrt( integrate(elements(mesh1), _Q<2*Order>(), (idv(u1)-idv(u1old))*(idv(u1)-idv(u1old)) ).evaluate()(0,0) );
	if (norm1!=0)  err1=err1/norm1;

	//U2old=U2;
	u2old=u2;
	p2old=p2;
	this->problem2(U2,U1);
	norm2 =math::sqrt( integrate( elements(mesh2), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );
	err2=math::sqrt( integrate(elements(mesh2), _Q<2*Order>(), (idv(u2)-idv(u2old))*(idv(u2)-idv(u2old)) ).evaluate()(0,0) );
	if (norm2!=0)  err2=err2/norm2;

	erreur=err1+err2;
	std::cout<<"\n erreur = "<<erreur<<" ; err1 = "<<err1<<" ; err2 = "<<err2;

	//exportResults( u1, u2, cpt);
	exportResults( U1, U2, cpt);

	++cpt;

      }

#endif
    //AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );

    //  u1= project( Xh1, elements(mesh1),  SolAnalytique );

    //u2= project( Xh2, elements(mesh2),  SolAnalytique );

    //exportResults( U1, U2, 0);

  }


  //--//////--//////--//////--//////-------------------------------------------------------------------------------------------------//
  //-//--//--//--//--//--//--//--//---------------------------------------------------------------------------------------//
  ////--//////--//////--//////--////------------------------------------------------------------------------------------------------//

  template<uint Dim,uint Order,uint OrderGeo>
  void
  Localisation<Dim,Order,OrderGeo>::problem1(element_u_type& u1 , element_l_type& l ) {
#if 0
    using namespace Life::vf;

    //----------------------------------------------------------------------------------------------------------//

    space_u_ptrtype Xh=u1.functionSpace();
    mesh_ptrtype mesh=Xh->mesh();

    element_u_type v(Xh,"v");

    sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );
    vector_ptrtype B( M_backend->newVector( Xh ) );
    A1->zero();
    B1->zero();

    //----------------------------------------------------------------------------------------------------------//

    value_type pi = M_PI;
    value_type nu= this->vm()["nu"].template as<value_type>();
    value_type alphaX =  this->vm()["alphaX"].template as<value_type>();
    value_type alphaY =  this->vm()["alphaY"].template as<value_type>();

    value_type penaldir = this->vm()["penaldir"].template as<double>();

    AUTO( Alpha, vec(cst(alphaX),cst(alphaY)) );
    AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
    AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

    //AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);
    AUTO( f,  -nu*lap_Uex);

    //----------------------------------------------------------------------------------------------------------//

    form2( Xh, Xh, A, _init=true ) =
      integrate( elements(mesh), _Q<2*Order>(),
		 nu*gradt(u1)*trans(grad(v)) );

    /*  form2( Xh, Xh, A) +=
	integrate( elements(mesh), _Q<2*Order>(),
	(gradt(u1)*Alpha)*id(v) );
    */

    form2( Xh, Xh, A ) +=
      integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		 id(v)*idt(p1)
		 +id(q)*idt(u1) );

    form2( Xh, Xh, A ) +=
      integrate( boundaryfaces(mesh), _Q<2*Order>(),
		 - nu*(gradt(u1)*N())*id(v));

    form2( Xh, Xh, A ) +=
      integrate( markedfaces(mesh,mesh->markerName("Bord")), _Q<2*Order>(),
		 //- nu*(grad(v)*N())*idt(u1)
		 + penaldir*id(v)*idt(u1)/hFace());

    form2( Xh, Xh, A ) +=
      integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		 //- nu*(grad(v)*N())*idt(u1)
		 + penaldir*id(q)*idt(p1)/hFace());


    A->close();


    form1( Xh, B, _init=true ) =
      integrate( elements(mesh), _Q<7>(),
		 f*id(v) );

    form1( Xh, B ) +=
      integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		 + penaldir*id(q)*idv(p2)/hFace()
		 -id(q)*idv(u2) );



#if 0
    form2( Xh, Xh, A ) +=
      //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
      on( markedfaces(mesh, mesh->markerName("Bord")) , U1, B, cst(0.) );

    //  form2( Xh, Xh, A ) +=
    //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
    //    on( markedfaces(mesh, mesh->markerName("Bord")) , p1, B, cst(0.) );
#endif

    this->solve(A,U1,B,false);

#endif
  }

  template<uint Dim,uint Order,uint OrderGeo>
  void
  Localisation<Dim,Order,OrderGeo>::problem2(element_u_type& u1 , element_l_type& l ) {
#if 0
    using namespace Life::vf;

    //----------------------------------------------------------------------------------------------------------//

    space_ptrtype Xh=U1.functionSpace();
    mesh_ptrtype mesh=Xh->mesh();

    element_type V(Xh,"v");
    element_0_type v = V.template element<0>();
    element_1_type q = V.template element<1>();

    element_0_type u1 = U1.template element<0>();
    element_1_type p1 = U1.template element<1>();

    element_0_type u2 = U2.template element<0>();
    element_1_type p2 = U2.template element<1>();

    sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );
    vector_ptrtype B( M_backend->newVector( Xh ) );
    A1->zero();
    B1->zero();

    //----------------------------------------------------------------------------------------------------------//

    value_type pi = M_PI;
    value_type nu= this->vm()["nu"].template as<value_type>();
    value_type alphaX =  this->vm()["alphaX"].template as<value_type>();
    value_type alphaY =  this->vm()["alphaY"].template as<value_type>();

    value_type penaldir = this->vm()["penaldir"].template as<double>();

    AUTO( Alpha, vec(cst(alphaX),cst(alphaY)) );
    AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
    AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

    //  AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);
    AUTO( f,  -nu*lap_Uex);

    //----------------------------------------------------------------------------------------------------------//
    
    form2( Xh, Xh, A, _init=true ) =
      integrate( elements(mesh), _Q<2*Order>(),
		 nu*gradt(u1)*trans(grad(v)) );

    /*  form2( Xh, Xh, A) +=
	integrate( elements(mesh), _Q<2*Order>(),
	(gradt(u1)*Alpha)*id(v) );
    */

    form2( Xh, Xh, A ) +=
      integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		 -id(v)*idt(p1) );

    form2( Xh, Xh, A ) +=
      integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		 -id(q)*idt(u1) );

    form2( Xh, Xh, A ) +=
      integrate( boundaryfaces(mesh), _Q<2*Order>(),
		 - nu*(gradt(u1)*N())*id(v));

    form2( Xh, Xh, A ) +=
      integrate( markedfaces(mesh,mesh->markerName("Bord")), _Q<2*Order>(),
		 //- nu*(grad(v)*N())*idt(u1)
		 + penaldir*id(v)*idt(u1)/hFace());

    form2( Xh, Xh, A ) +=
      integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		 + penaldir*id(q)*idt(p1)/hFace());

    //----------------------------------------------------------------------------------------------------------//

    form1( Xh, B, _init=true ) =
      integrate( elements(mesh), _Q<7>(),
		 f*id(v) );

    form1( Xh, B ) +=
      integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		 + penaldir*id(q)*idv(p2)/hFace()
		 +id(q)*idv(u2) );

    //----------------------------------------------------------------------------------------------------------//



    A->close();

#if 0
    form2( Xh, Xh, A ) +=
      //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
      on( markedfaces(mesh, mesh->markerName("Bord")) , U1, B, cst(0.) );

    //  form2( Xh, Xh, A ) +=
    //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
    //    on( markedfaces(mesh, mesh->markerName("Bord")) , p1, B, cst(0.) );
#endif

    this->solve(A,U1,B,false);
#endif
  }



  template<uint Dim,uint N,uint OrderGeo>
  void
  Localisation<Dim,N,OrderGeo>::solve( sparse_matrix_ptrtype const& D,
				       element_u_type& u,
				       vector_ptrtype const& F,
				       bool is_sym )
  {

    backend_ptrtype M_backend2(backend_type::build( this->vm() ));
  
    vector_ptrtype U( M_backend2->newVector( u.functionSpace() ) );
    M_backend2->solve( D, D, U, F, false );
    u = *U;
  } // Test_Schwarz::solve
 

  template<uint Dim, uint N, uint OrderGeo>
  template<typename IteratorRange, typename QuadratureT, typename ExprT>
  void
  Localisation<Dim,N,OrderGeo>::integrateLocal( IteratorRange const& itrange,
						QuadratureT const& im,
						ExprT const& expr,
						vector_ptrtype& rhs )
  {

#if 0
    mesh_ptrtype localmesh( new mesh_type );
    mesh_ptrtype globalmesh( new mesh_type );


    typedef space_u_type functionspace_type;
    space_u_ptrtype Uh = space_u_type::New( localmesh );
    typedef element_u_type element_type;  

    //
    // some typedefs
    //
    typedef typename mesh_type::element_type geoelement_type;
    typedef typename mesh_type::element_iterator mesh_element_iterator;
    typedef typename mesh_type::gm_type gm_type;
        
    // local
    typedef typename gm_type::template Context<ExprT::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, geoelement_type> local_gmc_type;
    typedef boost::shared_ptr<local_gmc_type> local_gmc_ptrtype;
        
    // global
    typedef typename gm_type::template Context<vm::POINT, geoelement_type> global_gmc_type;
    typedef boost::shared_ptr<global_gmc_type> global_gmc_ptrtype;
        
    //
    // Precompute some data in the reference element for
    // geometric mapping and reference finite element
    //
    typedef typename gm_type::precompute_type pc_type;
    typedef typename gm_type::precompute_ptrtype pc_ptrtype;
        
    //typedef typename QuadratureT::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex>::type quad_type;
        
    //    QuadMapped<quad_type> qm;
    //typedef typename QuadMapped<quad_type>::permutation_type permutation_type;
    //typename QuadMapped<quad_type>::permutation_points_type ppts( qm( im ) );
        
    QuadMapped<QuadratureT> qm;
    typedef typename QuadMapped<QuadratureT>::permutation_type permutation_type;
    typename QuadMapped<QuadratureT>::permutation_points_type ppts( qm( im ) );
        
    std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( im.nFaces() );
        
        
    //    typedef typename quad_type::face_quadrature_type face_im_type;
    typedef typename QuadratureT::face_quadrature_type face_im_type;
    //typedef typename QuadratureT::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex>::type::face_quadrature_type face_im_type;
        
    //__typeof__(im(__face_id_in_elt_0 ) ) im_face ( im(__face_id_in_elt_0 ) );
    std::vector<face_im_type> face_ims( im.nFaces() );
        
    for ( uint16_type __f = 0; __f < im.nFaces(); ++__f )
      {
	face_ims[__f] = im.face( __f );
                
	for( permutation_type __p( permutation_type::IDENTITY );
	     __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
	  {
	    __geopc[__f][__p] = pc_ptrtype(  new pc_type( localmesh->gm(), ppts[__f].find(__p)->second ) );
	  }
      }
        
    Log() << "[FBM] integrate Local\n";
        
    typedef typename boost::tuples::template element<1, IteratorRange>::type iterator_type;
    iterator_type local_it = itrange.template get<1>();
    iterator_type local_en = itrange.template get<2>();
        
    Log() << "[FBM::integrateLocal] n elts: " << std::distance( local_it, local_en ) << "\n";
    uint16_type __face_id_in_elt_0 = local_it->pos_first();
        
    // check that we have elements to iterate over
    if ( local_it == local_en )
      return;
        
    // get the geometric mapping associated with element 0
    local_gmc_ptrtype local_c0( new local_gmc_type( localmesh->gm(),
						    local_it->element( 0 ),
						    __geopc,//__geopc[__face_id_in_elt_0][local_it->element(0).permutation( __face_id_in_elt_0 )],
						    __face_id_in_elt_0 ) );
        
    uint16_type nptsface = im.nPointsOnFace();
        
    // locate the points in the global mesh
    typename mesh_type::Inverse meshinv( globalmesh );
        
    size_type nelts = std::distance(local_it,local_en);
    std::vector<double> expr_value( nelts*nptsface );
    for( size_type ide = 0; local_it != local_en; ++ local_it, ++ide )
      {
	//Log() << "[FBM] local element " << local_it->id() << "\n";
	// generate the points in the real space
	uint16_type __face_id_in_elt_0 = local_it->pos_first();
	local_c0->update( local_it->element(0),
			  //__geopc[__face_id_in_elt_0][local_it->element(0).permutation(__face_id_in_elt_0)],
			  __face_id_in_elt_0 );
                
	typedef fusion::map<fusion::pair<vf::detail::gmc<0>, local_gmc_ptrtype> > map_gmc_type;
                
	typedef typename ExprT::template tensor<map_gmc_type> eval_expr_type;
	typedef boost::shared_ptr<eval_expr_type> eval_expr_ptrtype;
	map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( local_c0 ) );
	eval_expr_ptrtype eval_expr( new eval_expr_type( expr, mapgmc ) );
                
                
	for( int q = 0; q < nptsface; ++ q )
	  {
	    size_type idq = nptsface*ide+q;
	    // evaluate expression and store it
	    expr_value[idq] = face_ims[__face_id_in_elt_0].weight(q)*eval_expr->evalq( 0, 0, q )*local_c0->J(q)*local_c0->normalNorm(q);

	    // insert quadrature point in inverse mesh data structure
	    meshinv.addPointWithId( local_c0->xReal( q ), idq, 0 );
	  }
      }
    meshinv.distribute();
        
        
        
    std::vector<bool> dof_done( nelts*nptsface );
    std::fill( dof_done.begin(), dof_done.end(), false );
    std::vector<boost::tuple<size_type,uint16_type > > itab;
    matrix_node<value_type>::type pts( mesh_type::nDim, 1 );
        
        
    mesh_element_iterator globalit = globalmesh->beginElementWithProcessId( Application::processId() );
    mesh_element_iterator globalen = globalmesh->endElementWithProcessId( Application::processId() );
        
    // basis
    typedef typename functionspace_type::basis_type basis_type;
    basis_type const* global_basis = Uh->basis().get();
    pc_ptrtype global_geopc( new pc_type( Uh->gm(), global_basis->dual().points() ) );
    global_gmc_ptrtype global_gmc( new global_gmc_type( Uh->gm(), *globalit, global_geopc ) );
    typedef typename basis_type::template Context<global_gmc_type::context,basis_type,gm_type,geoelement_type> basiscontext_type;
    typedef boost::shared_ptr<basiscontext_type> basiscontext_ptrtype;
        
        
    for( ; globalit != globalen; ++ globalit )
      {
	global_gmc->update( *globalit/*, global_geopc*/ );
	meshinv.pointsInConvex( globalit->id(), itab );
	if (itab.size() == 0)
	  continue;
                
	for (size_type i = 0; i < itab.size(); ++i)
	  {
	    // get dof id in target dof table
	    size_type q = boost::get<0>( itab[i] );
	    if ( !dof_done[q] )
	      {
		dof_done[q]=true;
		ublas::column( pts, 0 ) = meshinv.referenceCoords()[q];
                                
		global_geopc->update( pts );
		global_gmc_ptrtype global_gmccc( new global_gmc_type( Uh->gm(), *globalit, global_geopc ));
		//global_gmc->update( *globalit/*, global_geopc*/ );//Y a un souci?
		typedef typename element_type::pc_type elementpc_type;
		typedef typename element_type::pc_ptrtype elementpc_ptrtype;
		elementpc_ptrtype pc( new elementpc_type( Uh->fe(), pts ) );
		basiscontext_ptrtype global_bc( new basiscontext_type( Uh->fe(), global_gmccc, pc ) );
		for (uint16_type loc_ind=0; loc_ind < basis_type::nLocalDof ; loc_ind++)
		  {
		    for ( uint16_type comp = 0;comp < basis_type::nComponents;++comp )
		      {
			size_type globaldof = boost::get<0>(Uh->dof()->localToGlobal(globalit->id(), loc_ind, comp ));
                                                
			rhs->add( globaldof, expr_value[q]*global_bc->id( loc_ind, 0, 0, 0 ) );
		      }
		  }
                                
	      }
	  }
      } // global_it
        
        // make sure that all quad points have been visited
    for( size_type q = 0; q < dof_done.size(); ++q )
      {
	if ( dof_done[q] != true )
	  {
	    Log() << "quad point " << q << " has not been found in global mesh\n";
	  }
      }
    return;
#endif
  }


} // end namespace MortarMethod


/**
 * main function: entry point of the program
 */
int
main( int argc, char** argv )
{
  //dim, N, OrderGeo
  MortarMethod::Localisation<2,2,1> moncode( argc, argv,
					     MortarMethod::makeAbout(),
					     MortarMethod::makeOptions()
					     );

  moncode.runMultiDomain();

}



